// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtRadialResampleFilter.C - Test
// ************************************************************************* //
#include <float.h>

#include <avtRadialResampleFilter.h>
#include <avtParallel.h>
#include <avtCellLocatorBIH.h>

#include <DebugStream.h>

#include <vtkCellData.h>
#include <vtkCellType.h>
#include <vtkPoints.h>
#include <vtkPointData.h>
#include <vtkDoubleArray.h>
#include <avtWorldSpaceToImageSpaceTransform.h>

// ****************************************************************************
//  Method: avtRadialResampleFilter constructor
//
//  Programmer: griffin28 -- generated by xml2avt
//  Creation:   Tue May 20 13:15:11 PST 2014
//
// ****************************************************************************

avtRadialResampleFilter::avtRadialResampleFilter()
{
}


// ****************************************************************************
//  Method: avtRadialResampleFilter destructor
//
//  Programmer: griffin28 -- generated by xml2avt
//  Creation:   Tue May 20 13:15:11 PST 2014
//
//  Modifications:
//
// ****************************************************************************

avtRadialResampleFilter::~avtRadialResampleFilter()
{
}


// ****************************************************************************
//  Method:  avtRadialResampleFilter::Create
//
//  Programmer: griffin28 -- generated by xml2avt
//  Creation:   Tue May 20 13:15:11 PST 2014
//
// ****************************************************************************

avtFilter *
avtRadialResampleFilter::Create()
{
    return new avtRadialResampleFilter();
}


// ****************************************************************************
//  Method:      avtRadialResampleFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: Kevin Griffin
//  Creation:   Tue May 20 13:15:11 PST 2014
//
// ****************************************************************************

void
avtRadialResampleFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const RadialResampleAttributes*)a;
}


// ****************************************************************************
//  Method: avtRadialResampleFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtRadialResampleFilter with the given
//      parameters would result in an equivalent avtRadialResampleFilter.
//
//  Programmer: griffin28 -- generated by xml2avt
//  Creation:   Tue May 20 13:15:11 PST 2014
//
// ****************************************************************************
bool
avtRadialResampleFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(RadialResampleAttributes*)a);
}


// ****************************************************************************
//  Method: avtRadialResampleFilter::Execute
//
//  Purpose:
//      Sends the specified input and output through the RadialResample filter.
//
//  Returns: The output dataset.
//
//  Programmer: Kevin Griffin
//  Creation:   Tue May 20 13:15:11 PST 2014
//
//  Modifications:
//    Kevin Griffin, Fri Mar 24 18:06:06 PDT 2017
//    Fixed a memory leak and added logic to support point or cell data.
//
// ****************************************************************************
void
avtRadialResampleFilter::Execute()
{
    avtDataset_p output = GetTypedOutput();
    double bounds[6] = { 0, 0, 0, 0, 0, 0 };
    GetBounds(bounds);

    debug4 << "Resampling over space: " << bounds[0] << ", " << bounds[1]
           << ": " << bounds[2] << ", " << bounds[3] << ": " << bounds[4]
           << ", " << bounds[5] << endl;
    
    //
    // Our resampling leaves some invalid values in the data range.  The
    // easiest way to bypass this is to get the data range from the input and
    // pass it along (since resampling does not change it in theory).
    //
    double range[2];
    if (GetInput()->GetInfo().GetAttributes().ValidActiveVariable())
    {
        GetDataExtents(range);
        output->GetInfo().GetAttributes().GetDesiredDataExtents()->Set(range);
    }

    //
    // If there are no variables, then just create the mesh and exit.
    //
    bool thereAreNoVariables = (GetInput()->GetInfo().GetAttributes().GetNumberOfVariables() <= 0);
    
    if(thereAreNoVariables)
    {
        if(PAR_Rank() == 0)
        {
            vtkDataSet *rrg;    /* Radial Resample Grid */
            if(atts.GetIs3D())
            {
                rrg = Create3DRadialGrid(bounds,
                                        atts.GetCenter(),
                                        atts.GetMinAzimuth(),
                                        atts.GetMaxAzimuth(),
                                        atts.GetDeltaAzimuth(),
                                        atts.GetMinTheta(),
                                        atts.GetMaxTheta(),
                                        atts.GetDeltaTheta(),
                                        atts.GetRadius(),
                                        atts.GetDeltaRadius(),
                                        atts.GetIsFast());
            } else
            {
                rrg = Create2DRadialGrid(bounds,
                                        atts.GetCenter(),
                                        atts.GetMinTheta(),
                                        atts.GetMaxTheta(),
                                        atts.GetDeltaTheta(),
                                        atts.GetRadius(),
                                        atts.GetDeltaRadius(),
                                        atts.GetIsFast());
            }
        
        
            avtDataTree_p tree = new avtDataTree(rrg, 0);
            rrg->Delete();
            SetOutputDataTree(tree);
        }
        else
        {
            //
            // Putting in a NULL data tree can lead to seg faults, etc.
            //
            avtDataTree_p dummy = new avtDataTree();
            SetOutputDataTree(dummy);
        }
        
        return;
    }
    else // Have variables
    {
        vtkDataSet *rrg;
        
        if(atts.GetIs3D())
        {
            rrg = Create3DRadialGrid(bounds,
                                    atts.GetCenter(),
                                    atts.GetMinAzimuth(),
                                    atts.GetMaxAzimuth(),
                                    atts.GetDeltaAzimuth(),
                                    atts.GetMinTheta(),
                                    atts.GetMaxTheta(),
                                    atts.GetDeltaTheta(),
                                    atts.GetRadius(),
                                    atts.GetDeltaRadius(),
                                    atts.GetIsFast());
        } else
        {
            rrg = Create2DRadialGrid(bounds,
                                    atts.GetCenter(),
                                    atts.GetMinTheta(),
                                    atts.GetMaxTheta(),
                                    atts.GetDeltaTheta(),
                                    atts.GetRadius(),
                                    atts.GetDeltaRadius(),
                                    atts.GetIsFast());
        }
        
        avtDataTree_p inputTree = GetInputDataTree();   /* get input data tree to obtain datasets */
        int nsets;                                      /* number of datasets */
        int i,j,k;                                      /* for loop indexes */
        int rrgNumOfPoints = rrg->GetNumberOfPoints();
        
        vtkDataSet **dataSets = inputTree->GetAllLeaves(nsets);
        debug5 << "nsets = " << nsets << endl;
        
        avtCellLocatorBIH **cellLocators = NULL;
        if(nsets > 0)
        {
            cellLocators = new avtCellLocatorBIH*[nsets];
            
            // Create cell locator for each dataset
            for(k=0; k<nsets; k++)
            {
                cellLocators[k] = new avtCellLocatorBIH(dataSets[k]);
            }
        }
                
        int numOfVars = GetInput()->GetInfo().GetAttributes().GetNumberOfVariables();
        vtkDoubleArray **vars = new vtkDoubleArray*[numOfVars];     /* Resampled variable values generated by interpolating the cell values
                                                                        of the original dataset */
        std::vector<vtkDataObject::FieldAssociations> fieldAssocs(numOfVars);
        avtInterpolationWeights lastWeights;                        /* Interpolation weights */
        vtkIdType lastCellId;                                       /* Id of the original cell containing the resampled point */
        double pos[3];                                              /* Resampled point */
        
        for(i=0; i<numOfVars; i++)
        {
            std::string varName = GetInput()->GetInfo().GetAttributes().GetVariableName(i);
            
            int numVarDataComps = GetInput()->GetInfo().GetAttributes().GetVariableDimension(varName.c_str());
            debug5 << "varName: " << varName << " Number of Components: " << numVarDataComps <<  endl;
            
            // Get the scalar values from the current grid (inputArray) and
            // interpolate for the new grid (var)
            vars[i] = vtkDoubleArray::New();
            vars[i]->SetName(varName.c_str());
            vars[i]->SetNumberOfTuples(rrgNumOfPoints);
            vars[i]->SetNumberOfComponents(numVarDataComps);
            
            /* orginal variable value */
            double *inputArrayDataVals = new double[numVarDataComps];
            /* variable values after weight applied */
            double *outputArrayDataVals = new double[numVarDataComps];
            
            for(j=0; j<rrgNumOfPoints; j++)
            {
                rrg->GetPoint(j, pos);
               
                for(int compIdx=0; compIdx<numVarDataComps; compIdx++)
                {
                    outputArrayDataVals[compIdx] = 0;
                }
                
                for(k=0; k<nsets; k++)  // loop over all the datasets (leaves)
                {
                    vtkDataSet *dataSet = dataSets[k];
                    vtkDataArray *inputArray = dataSet->GetPointData()->GetArray(varName.c_str());
                    if(inputArray != NULL)
                    {
                        fieldAssocs[i] = vtkDataObject::FIELD_ASSOCIATION_POINTS;
                    }
                    else
                    {
                        inputArray = dataSet->GetCellData()->GetArray(varName.c_str());
                        fieldAssocs[i] = vtkDataObject::FIELD_ASSOCIATION_CELLS;
                    }
                    
                    if(inputArray != NULL)
                    {
                        lastCellId = cellLocators[k]->FindCell(pos, &lastWeights, true);
                                                
                        if(lastCellId != -1)
                        {
                            for (avtInterpolationWeights::const_iterator wi=lastWeights.begin(); wi!=lastWeights.end(); ++wi)
                            {
                                inputArray->GetTuple(wi->i, inputArrayDataVals);
                    
                                for(int compIdx=0; compIdx < numVarDataComps; compIdx++)
                                {
                                    // Multiply each component of the old scalar by the weight and update the new scalar value
                                    outputArrayDataVals[compIdx] += wi->w * inputArrayDataVals[compIdx];
                                }
                            }
                            
                            break;
                        }
                    }
                }
                
                vars[i]->SetTuple(j, outputArrayDataVals);
            }
            
            // Reduction if running in parallel: MPI_Reduce
            int count = numVarDataComps * rrgNumOfPoints;
            double *outResult = (double *) vars[i]->GetVoidPointer(0);
            Collect(outResult, count);
            
            // Memory Cleanup
            delete [] inputArrayDataVals;
            delete [] outputArrayDataVals;
        }
        
        if(PAR_Rank() == 0)
        {
            // Attach variables to our radial grid
            for(i=0; i<numOfVars; i++)
            {
                if(fieldAssocs[i] == vtkDataObject::FIELD_ASSOCIATION_POINTS)
                {
                    rrg->GetPointData()->AddArray(vars[i]);
                }
                else
                {
                    rrg->GetCellData()->AddArray(vars[i]);
                }
                
                if(i == 0)  // Assuming the first variable is the primary variable
                {
                    if(vars[i]->GetNumberOfComponents() == 1)
                    {
                        if(fieldAssocs[i] == vtkDataObject::FIELD_ASSOCIATION_POINTS)
                        {
                            rrg->GetPointData()->SetActiveScalars(vars[i]->GetName());
                        }
                        else
                        {
                            rrg->GetCellData()->SetActiveScalars(vars[i]->GetName());
                        }
                    }
                    else if(vars[i]->GetNumberOfComponents() == 3)
                    {
                        if(fieldAssocs[i] == vtkDataObject::FIELD_ASSOCIATION_POINTS)
                        {
                            rrg->GetPointData()->SetActiveVectors(vars[i]->GetName());
                        }
                        else
                        {
                            rrg->GetCellData()->SetActiveVectors(vars[i]->GetName());
                        }
                    }
                }
            }
            
            avtDataTree_p tree = new avtDataTree(rrg, 0);
            rrg->Delete();            
            SetOutputDataTree(tree);
        }
        else
        {
            //
            // Putting in a NULL data tree can lead to seg faults, etc.
            //
            avtDataTree_p dummy = new avtDataTree();
            SetOutputDataTree(dummy);
        }
        
        // Memory Cleanup
        for(i=0; i<numOfVars; i++)
        {
            vars[i]->Delete();
        }
        
        delete [] vars;
        
        if(cellLocators != NULL)
        {
            for(k=0; k<nsets; k++)
            {
                delete cellLocators[k];
            }
        
            delete [] cellLocators;
        }
    }
}

// ****************************************************************************
//  Method: avtResampleFilter::GetBounds
//
//  Purpose:
//      Obtains the bounds of the resampled volume.  This could come from
//      attributes, or the existing spatial extents of the input.
//
//  Arguments:
//      bounds       Output array.  Format is min/max X, then m/m Y, m/m Z.
//
//  Returns: whether or not these specify 3-dimensional bounds.
//
//  Programmer: Tom Fogal
//  Creation:   June 23, 2009
//
//  Modifications:
//
//    Hank Childs, Thu Aug 26 13:47:30 PDT 2010
//    Change extents names.
//
//    Hank Childs, Tue Nov 30 21:54:43 PST 2010
//    Remove const qualification.
//
//    Kevin Griffin, Tue May 13 15:30:12 PDT 2014
//    Modified for use in this class by removing the
//    atts.GetUseBounds() conditional
//
// ****************************************************************************
bool avtRadialResampleFilter::GetBounds(double bounds[6])
{
    bool is3D = true;
    const avtDataAttributes &datts = GetInput()->GetInfo().GetAttributes();
    avtExtents *exts = datts.GetDesiredSpatialExtents();
    
    if (exts->HasExtents())
    {
        exts->CopyTo(bounds);
    }
    else
    {
        GetSpatialExtents(bounds);
    }

    if (fabs(bounds[4]) < 1e-100 && fabs(bounds[5]) < 1e-100)
    {
        is3D = false;
        bounds[5] += 0.1;
    }
    return is3D;
}

// ****************************************************************************
//  Method: avtRadialResampleFilter::CreateFast2DRadialGrid
//
//  Purpose:
//      Creates a circular grid using a structured grid for efficiency. The
//      grid will have a small hole (0.001 radius) since structured grids can
//      only use cells that are quadrilaterals. 
//
//  Arguments:
//      dims        The x and y dimensions of the structured grid (z = 1).
//      center      The center of the grid
//      startTheta  The start angle.
//      deltaTheta  The angle delta.
//      deltaRadius The radius delta.
//
//  Programmer:   Kevin Griffin
//  Creation:     May 13, 2014
//
//  Modifications:
//
// ****************************************************************************
vtkStructuredGrid *
avtRadialResampleFilter::CreateFast2DRadialGrid(int dims[3], const float *center, float startTheta, float deltaTheta, float deltaRadius)
{
    int i, j;
    int jOffset, offset;
    float theta, radius;
    float coords[3];
    coords[2] = center[2];
    
    vtkStructuredGrid *sgrid = vtkStructuredGrid::New();
    sgrid->SetDimensions(dims);
    
    vtkPoints *points = vtkPoints::New();
    points->Allocate(dims[0]*dims[1]);
    
    for(j=0; j<dims[1]; j++)
    {
        jOffset = j * dims[0];
        for(i=0; i<dims[0]; i++)
        {
            if(j != 0)
            {
                radius = j*deltaRadius;
            } else
            {
                radius = 0.001;
            }
            
            theta = (startTheta + (i * deltaTheta)) * M_PI/180;   // Convert to radians
            coords[0] = center[0] + radius * cos(theta);
            coords[1] = center[1] + radius * sin(theta);
            
            offset = i + jOffset;
            points->InsertPoint(offset, coords);
        }
    }
    
    sgrid->SetPoints(points);
    sgrid->Squeeze();
    points->Delete();
    
    return sgrid;
}

// ****************************************************************************
//  Method: avtRadialResampleFilter::CreateNormal2DRadialGrid
//
//  Purpose:
//      Creates a circular grid using an unstructured grid with a combination
//      of triangles and quadrilaterals. 
//
//  Arguments:
//      dims        The x and y dimensions of the structured grid (z = 1).
//      center      The center of the grid
//      startTheta  The start angle.
//      deltaTheta  The angle delta.
//      deltaRadius The radius delta.
//
//  Programmer:   Kevin Griffin
//  Creation:     May 13, 2014
//
//  Modifications:
//
// ****************************************************************************
vtkUnstructuredGrid *
avtRadialResampleFilter::CreateNormal2DRadialGrid(int dims[3], const float *center, float startTheta, float deltaTheta, float deltaRadius)
{
    int i, j;
    int jOffset, baseOffset, offset;
    float theta, radius;
    float coords[] = {0.0, 0.0, 0.0};
    
    vtkPoints *points = vtkPoints::New();
    points->Allocate(dims[0] * dims[1]);
    
    // Start with center point
    coords[0] = center[0];
    coords[1] = center[1];
    coords[2] = center[2];
    points->InsertPoint(0, coords);
    
    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
    ugrid->Allocate();
    
    for(j=1; j<dims[1]; j++)
    {
        jOffset = j != 1 ? (j-1) * dims[0] + 1 : 1;
        for(i=0; i<dims[0]; i++)
        {
            radius = j * deltaRadius;
            offset = i + jOffset;
            
            // Calculate Point
            theta = (startTheta + (i * deltaTheta)) * M_PI/180;   // Convert to radians
            coords[0] = center[0] + radius * cos(theta);
            coords[1] = center[1] + radius * sin(theta);
            points->InsertPoint(offset, coords);
            
            if(i != 0)
            {
                if(j != 1)
                {    // Quadrilateral Cells
                    baseOffset = i + ((j-2) * dims[0] + 1);
                    vtkIdType pts[] = {baseOffset, baseOffset-1, offset-1, offset};
                    ugrid->InsertNextCell(VTK_QUAD, 4, pts);
                }
                else
                {    // Triangle Cells
                    vtkIdType pts[] = {0, offset-1, offset};
                    ugrid->InsertNextCell(VTK_TRIANGLE, 3, pts);
                }
            } 
        }
    }
    
    ugrid->SetPoints(points);
    ugrid->Squeeze();
    
    points->Delete();
    
    return ugrid;
}

// ****************************************************************************
//  Method: avtRadialResampleFilter::Create2DRadialGrid
//
//  Purpose:
//      Creates a circular grid parameterized by the start and stop angle.
//      If isFast is true a more efficient structured grid will be used
//      leaving a small hole in the grid. This is due to structured
//      grids only being able to compose cells that are either quadrilaterals
//      (2D) or hexahedrons (3D).Otherwise an unstructured grid will be used
//      leaving no hole in the grid.
//
//  Arguments:
//      bounds      The bounds of the original dataset.
//      minTheta    The start angle.
//      maxTheta    The stop angle.
//      deltaTheta  The angle increment amount when going from minTheta
//                  to maxTheta
//      radius      The radius of the grid
//      deltaRadius The radius increment
//      isFast      If true use a structured grid, otherwise use an
//                  unstructured grid.
//
//  Programmer:   Kevin Griffin
//  Creation:     May 13, 2014
//
//  Modifications:
//
// ****************************************************************************
vtkDataSet *
avtRadialResampleFilter::Create2DRadialGrid(const double *bounds, const float *center, float minTheta, float maxTheta,
                                            float deltaTheta, float radius, float deltaRadius, bool isFast)
{
    int dims[3];
    dims[0] = static_cast<int>(std::fabs((maxTheta - minTheta)/deltaTheta)) + 1;
    dims[1] = static_cast<int>(radius/deltaRadius) + 1;
    dims[2] = 1;
    
    if(isFast)
    {
        return CreateFast2DRadialGrid(dims, center, minTheta, deltaTheta, deltaRadius);
    }
    else
    {
        return CreateNormal2DRadialGrid(dims, center, minTheta, deltaTheta, deltaRadius);
    }
}

// ****************************************************************************
//  Method: avtRadialResampleFilter::CreateNormal3DRadialGrid
//
//  Purpose:
//      Creates a spherical grid parameterized by the azimuth, elevation, and
//      radius.
//
//  Arguments:
//      dims            The dimensions of the Spherical Grid
//      center          The center coordinates to use.
//      startAzimuth    The start angle in the XZ plane.
//      startElevation  The start angle in the XY plane.
//      deltaAzimuth    The angle step to use when traversing the start azimuth
//                      to the stop azimuth
//      delaElevation   The angle step to use when traversing from start elevation
//                      to stop elevation.
//      deltaRadius     The radius delta.
//
//  Programmer:   Kevin Griffin
//  Creation:     May 16, 2014
//
//  Modifications:
//
// ****************************************************************************
vtkUnstructuredGrid *
avtRadialResampleFilter::CreateNormal3DRadialGrid(int dims[3], const float *center, float startAzimuth, float startElevation,
                                                  float deltaAzimuth, float deltaElevation, float deltaRadius)
{
    int i, j, k;
    int kOffset, jOffset;
    int offset, baseOffset;
    int kOffset1, jOffset1;
    int offset1, baseOffset1;
    float theta, phi, radius;
    float coords[3];
    
    vtkPoints *points = vtkPoints::New();
    points->Allocate(dims[0]*dims[1]*dims[2]);
    
    // Add center point
    coords[0] = center[0];
    coords[1] = center[1];
    coords[2] = center[2];
    points->InsertPoint(0, coords);
    
    vtkUnstructuredGrid *ugrid = vtkUnstructuredGrid::New();
    ugrid->Allocate();
    
    for(k=0; k<dims[2]; k++)
    {
        kOffset = k * dims[0] * dims[1];
        kOffset1 = (k-1) * dims[0] * dims[1];
        for(j=1; j<dims[1]; j++)
        {
            jOffset = j != 1 ? (j-1) * dims[0] + 1 : 1;
            jOffset1 = (j-2) * dims[0] + 1;
            for(i=0; i<dims[0]; i++)
            {
                radius = j*deltaRadius;
                offset = i + jOffset + kOffset;
                
                phi = (startAzimuth + (k * deltaAzimuth)) * M_PI/180;
                theta = (startElevation + (i * deltaElevation)) * M_PI/180;
                
                coords[0] = (center[0] + radius * cos(theta) * cos(phi));
                coords[1] = center[1] + radius * sin(theta);
                coords[2] = center[2] + radius * cos(theta) * sin(phi);
                points->InsertPoint(offset, coords);
                
                if(k != 0)
                {
                    if(i != 0)
                    {
                        offset1 = i + jOffset + kOffset1;
                        
                        if(j != 1)  // Hexahedron
                        {
                            baseOffset = i + jOffset1 + kOffset;
                            baseOffset1 = i + jOffset1 + kOffset1;
                            
                            vtkIdType pts[] = {baseOffset1-1, baseOffset1, offset1, offset1-1, baseOffset-1, baseOffset, offset, offset-1};
                            ugrid->InsertNextCell(VTK_HEXAHEDRON, 8, pts);
                        }
                        else        // Pyramid
                        {
                            vtkIdType pts[] = {offset1-1, offset1, offset, offset-1, 0};
                            ugrid->InsertNextCell(VTK_PYRAMID, 5, pts);
                        }
                    }
                }
            }
        }
    }
    
    ugrid->SetPoints(points);
    ugrid->Squeeze();
    points->Delete();
    
    return ugrid;
}

// ****************************************************************************
//  Method: avtRadialResampleFilter::CreateFast3DRadialGrid
//
//  Purpose:
//      Creates a spherical grid parameterized by the azimuth, elevation, and
//      radius. A more efficient structured grid is used, however, a small hole
//      is left in the grid due to the fact that a structured grid can only
//      compose cells that are hexahedron. Use this method ONLY for a quick
//      estimation.
//
//  Arguments:
//      dims            The dimensions of the Spherical Grid
//      center          The center coordinates to use.
//      startAzimuth    The start angle in the XZ plane.
//      startElevation  The start angle in the XY plane.
//      deltaAzimuth    The angle step to use when traversing the start azimuth
//                      to the stop azimuth
//      delaElevation   The angle step to use when traversing from start elevation
//                      to stop elevation.
//      deltaRadius     The radius delta.
//
//  Programmer:   Kevin Griffin
//  Creation:     May 16, 2014
//
//  Modifications:
//
// ****************************************************************************
vtkStructuredGrid *
avtRadialResampleFilter::CreateFast3DRadialGrid(int dims[3], const float *center, float startAzimuth, float startElevation,
                                                float deltaAzimuth, float deltaElevation, float deltaRadius)
{
    int i, j, k;
    int kOffset;
    int jOffset;
    int offset;
    float theta;
    float phi;
    float radius;
    float coords[3];
    
    vtkStructuredGrid *sgrid = vtkStructuredGrid::New();
    sgrid->SetDimensions(dims);
    
    vtkPoints *points = vtkPoints::New();
    points->Allocate(dims[0]*dims[1]*dims[2]);
    
    for(k=0; k<dims[2]; k++)
    {
        kOffset = k * dims[0] * dims[1];
        for(j=0; j<dims[1]; j++)
        {
            jOffset = j * dims[0];
            for(i=0; i<dims[0]; i++)
            {
                if(j != 0)
                {
                    radius = j*deltaRadius;
                }
                else
                {
                    radius = 0.001;
                }
                
                phi = (startAzimuth + (k * deltaAzimuth)) * M_PI/180;
                theta = (startElevation + (i * deltaElevation)) * M_PI/180;
                
                coords[0] = center[0] + radius * cos(theta) * cos(phi);
                coords[1] = center[1] + radius * sin(theta);
                coords[2] = center[2] + radius * cos(theta) * sin(phi);
                
                offset = i + jOffset + kOffset;
                points->InsertPoint(offset, coords);
            }
        }
    }
    
    sgrid->SetPoints(points);
    points->Delete();
    
    return sgrid;
}

// ****************************************************************************
//  Method: avtRadialResampleFilter::Create3DRadialGrid
//
//  Purpose:
//      Creates a spherical grid parameterized by the azimuth (phi), elevation (theta),
//      and radius. If isFast is true a more efficient structured grid will be used
//      leaving a small hole in the center of the sphere. This is due to structured
//      grids only being able to compose cells that are either quadrilaterals (2D) or
//      hexahedrons (3D). Otherwise an unstructured grid will be used leaving no hole
//      to be filled.
//
//  Arguments:
//      bounds          The bounds of the original dataset.
//      center          The center coordinates to use.
//      minAzimuth      The start angle in the XZ plane.
//      maxAzimuth      The stop angle in the XZ plane.
//      deltaAzimuth    The angle step to use when traversing the start azimuth
//                      to the stop azimuth
//      minElevation    The start angle in the XY plane.
//      maxElevation    The stop angle in the XY plane
//      delaElevation   The angle step to use when traversing from minElevation
//                      to maxElevation.
//      radius          The radius of the sphere.
//      deltaRadius     The radius delta.
//      isFast          If true use a structured grid, otherwise use an
//                      unstructured grid.
//
//  Programmer:   Kevin Griffin
//  Creation:     May 16, 2014
//
//  Modifications:
//
// ****************************************************************************
vtkDataSet *
avtRadialResampleFilter::Create3DRadialGrid(const double *bounds, const float *center, float minAzimuth, float maxAzimuth, float deltaAzimuth,
                                            float minElevation, float maxElevation, float deltaElevation, float radius, float deltaRadius, bool isFast)
{
    int dims[3];    /* Dimensions of the Resampled Grid */
    
    dims[0] = static_cast<int>(std::fabs((maxElevation - minElevation)/deltaElevation)) + 1;
    dims[1] = static_cast<int>(radius/deltaRadius) + 1;
    dims[2] = static_cast<int>(std::fabs((maxAzimuth - minAzimuth)/deltaAzimuth)) + 1;
    
    if(isFast)
    {
        return CreateFast3DRadialGrid(dims, center, minAzimuth, minElevation, deltaAzimuth, deltaElevation, deltaRadius);
    }
    else
    {
        return CreateNormal3DRadialGrid(dims, center, minAzimuth, minElevation, deltaAzimuth, deltaElevation, deltaRadius);
    }
}

avtContract_p
avtRadialResampleFilter::ModifyContract(avtContract_p in_contract)
{
    avtContract_p rv = new avtContract(in_contract);
    resampleVarName = std::string(rv->GetDataRequest()->GetVariable());
    
    if (in_contract->GetDataRequest()->MayRequireZones() ||
        in_contract->GetDataRequest()->MayRequireNodes())
    {
        avtDataAttributes &data = GetInput()->GetInfo().GetAttributes();
        
        if (data.ValidActiveVariable())
        {
            if (data.GetCentering() == AVT_NODECENT)
            {
                rv->GetDataRequest()->TurnNodeNumbersOn();
            }
            else if (data.GetCentering() == AVT_ZONECENT)
            {
                rv->GetDataRequest()->TurnZoneNumbersOn();
            }
        }
        else
        {
            // canot determine variable centering, so turn on both
            // node numbers and zone numbers.
            rv->GetDataRequest()->TurnNodeNumbersOn();
            rv->GetDataRequest()->TurnZoneNumbersOn();
        }
    }
    
    return rv;
}

void
avtRadialResampleFilter::UpdateDataObjectInfo(void)
{
//    avtDataAttributes &inAtts   = GetInput()->GetInfo().GetAttributes();
    avtDataAttributes &outAtts = GetOutput()->GetInfo().GetAttributes();
    avtDataValidity   &outValidity = GetOutput()->GetInfo().GetValidity();
//    int spatialDim = inAtts.GetSpatialDimension();
//    
//    outAtts.SetTopologicalDimension(spatialDim);
    
    outValidity.InvalidateZones();
//    outValidity.InvalidateNodes();
    
    if(!resampleVarName.empty() && outAtts.ValidActiveVariable())
    {
        if(outAtts.GetVariableName() != resampleVarName)
        {
            outAtts.SetActiveVariable(resampleVarName.c_str());
        }
    }
    
    char params[200];
    snprintf(params, 200, "Min Elevation=%f Max Elevation=%f Delta=%f Radius=%f", atts.GetMinTheta(), atts.GetMaxTheta(), atts.GetDeltaTheta(), atts.GetRadius());
    outAtts.AddFilterMetaData("Radial Resample", params);
}
